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Abstract. We present an alternative method to estimate the numerical viscosity 
in simulations of astrophysical objects, which is based in the damping of 
fluid oscillations. We apply the method to general relativistic hydrodynamic 
simulations using spherical coordinates. We perform ID-spherical and 2D- 
axisymmetric simulations of radial oscillations in spherical systems. We calibrate 
flrst the method with simulations with physical bulk viscosity and study the 
differences between several numerical schemes. We apply the method to radial 
oscillations of neutron stars and we conclude that the main source of numerical 
viscosity in this case is the surface of the star. We expect that this method could 
be useful to compute the resolution requirements and limitations of the numerical 
simulations in different astrophysical scenarios in the future. 
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1. Introduction 

There is a number of scenarios in theoretical astrophysics which need to be modeUed 
by means of numerical simulations. The complexity found in these scenarios does not 
allow for a simple analytical study to explain the observations. In a subset of these 
scenarios the effects of general relativity are important and a Newtonian description of 
the system is not sufficient. Some examples of stellar size scenarios of this kind are: i) 
the collapse of massive stellar cores to form a neutron star, which eventually can lead 
to a supernova explosion, ii) for more massive stars, the collapse to a black hole and 
possibly an accretion disk, which can power a long-duration gamma-ray burst (GRB) 
according to the collapsar scenario (Woosley 1993) and could be related to the subclass 
of broad-lined type Ic supernovae (Modjaz et al 2008), and iii) the cooling of a rapidly 
rotating hot proto-neutron star (PNS) to form a magnetar (Thomson & Duncan 
1993). All three previous examples have in common that are quasi-spherical scenarios. 
Therefore it is natural to model them in numerical simulations using spherical polar 
coordinates. 

Those scenarios have some common points. The three scenarios include multi- 
scale physics in which a broad range of length and time scales play an important 
role. In the supernova core collapse case the length scales to cover range from 
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the size of the iron core (~ 1000 km) to the size of the PNS (~ 10 km). The 
smaller length-scale process dominating the dynamics ranges from ^ 1 km in the 
case of convectively unstable PNS formed after the collapse of the most common slow 
rotating progenitors (Miiller & Janka 1997), to 1 cm - 1 m if the magneto-rotational 
instability develops in the collapse of fast rotating progenitors (Cerda-Duran et al 
2007). Furthermore the mechanism for the explosion requires an appropriate modeling 
of the neutrino transport. The energy deposited by neutrinos behind the shock and 
the standing accretion-shock instability (SASI) probably arc crucial for a successful 
explosion (Marek & Janka 2009) . In the coUapsar scenario strong rotation is necessary 
for the formation of a GRB. It allows the generation of an accretion disk around the 
formed black hole with a low-density funnel along the rotation axis. The relativistic jet 
which is responsible for the GRB is powered either by MHD processes or by the energy 
deposition of annihilating neutrinos around the axis (see e.g. MacFadyen et al 2001). 
In the first case, magneto-rotational instabilities and dynamo processes arc most 
probably responsible for the amplification of the magnetic field during the collapse. 
The case of dynamo processes in rapidly rotating PNS resembles the modeling of the 
sun. The cooling timc-scalc is of the order of seconds while the rotation period which 
drives the dynamo processes can be as low as 1 ms reaching magnetic Reynolds number 
of ~ 10^^ (Thompson & Duncan 1993). The numerical modeling of all scenarios 
described above require very high resolution to resolve the MRI and the turbulence 
amplifying the magnetic field, full three-dimensional simulations to correctly capture 
the growth and saturation of the different types of instabilities (MRI, SASI, global 
MHD instabilities), and a time evolution much larger than the characteristic dynamical 
timescale of the system. Therefore the numerical methods intended to solve this 
problem should be as less dissipative as possible in order to guarantee that the 
numerical dissipation is smaller than the expected physical one or, in cases where 
the computational requirements do not allow this, the global properties of the system 
show convergence with increasing grid resolution. The numerical dissipation could be 
a serious limitation for the successful modeling of objects mentioned above. It is hence 
important to have the appropriate scaling with the resolution, which is a necessary 
property to have scalable numerical codes which can run efficiently in massive parallel 
computers. 

Most of the numerical codes performing simulations of these scenarios use Eulerian 
grids, explicit numerical schemes, and a mesh adapted to the problem, mostly with 
grids in spherical polar coordinates. Eulerian grid-based codes are better suited for 
these scenarios than Lagrangian methods (e.g. smoothed particle hydrodynamics) 
because they allow to use finite-volume conservative schemes. Eulerian methods allow 
for the correct treatment of arbitrarily high discontinuities and shocks in general 
relativity (Ibafiez et al 2000, Dimmelmeier et al 2002, Duez et al 2003, Shibata 2007, 
Baiotti et al 2007), even with magnetic fields (Gammie et al 2003, Komissarov et al 
2005, Anninos et al 2005, Anton et al 2006). Explicit methods are better suited for 
multidimensional simulations since they are computationally less expensive and easier 
to parallelize, although they have time-step limitations given by the CFL condition 
(Courant et al 1928, 1967). Spherical polar coordinates have several properties which 
make them appropriate to model the objects described above: i) are well adapted 
for quasi-spherical objects, ii) allow for accurate conservation of angular momentum, 
which is not true, in general, for Eulerian grids (sec e.g. Zink et al 2008 and Fragile 
et al 2009), iii) axisymmetry and spherical symmetry can be easily enforced, and iv) 
large radial domains can be covered using non equally spaced radial grids. These 
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numerical methods have been successfuUy apphed in ID (spherical symmetry) and 
2D (axisymmetry) simulations. However the extension to 3D in spherical coordinates 
suffers from severe time-step restrictions which render simulations unaffordable unless 
a special treatment of the central and polar regions is used (see. e.g Mller & Janka 
1997) 

Numerical dissipation effects are present in simulations using Eulerian grids due 
to the discretization of the equations that are solved. There are several methods 
to quantify the amount of numerical dissipation of a code based on the measure 
of the energy losses of the system. This has been standard practice since the first 
studies of hydrodynamic turbulence (e.g. Herring et al 1974) due to the necessity of 
resolving the physical dissipation scales. There are also studies of the decay of waves in 
hydrodynamics (Porter et al 1994) and MHD simulations (e.g. Simon & Hawley 2009). 
More recently the numerical dissipation has been estimated measuring the angular 
momentum transport by MHD turbulence (Fromang & Papaloizou 2007, Simon et al 
2009). All these methods allow for a simplified numerical setup in where the local 
dissipation properties of the numerical algorithms can accurately be estimated. We 
propose an alternative approach to measure the numerical dissipation that is suitable 
for global simulations of relativistic stars close to the equilibrium. 

The aim of this paper is to study the effects of numerical viscosity in simulations 
with spherical coordinates and study the influence of the grid resolution and the 
numerical scheme. In section [2] we describe the hydrodynamics equations including 
physical bulk viscosity in general relativity (GR), in section [3] we present a method to 
estimate numerical dissipation effects in a simplified test case. We apply this method 
in section |4] to estimate the numerical viscosity of oscillating neutron stars. We finish 
the paper in section [5] discussing the implications of our numerical results. If is not 
explicitly mentioned, we use units in which c = G = 1. Greek indices run from to 3 
and Latin indices from 1 to 3. 

2. GR hydrodynamics with bulk viscosity 

We use the 3+1 decomposition of the spacetime in which the metric reads 

ds^ = g^.^dx^'dx" = -c? dt^ + 7y (dx* + 13' dt){dx^ + P'^ dt), (1) 

where a, /3* and 7'-' are the lapse function, the shift vector and the spatial 3- 
metric respectively. In addition we consider the conformally flat condition (CFC) 
approximation (Isenberg 2008, Wilson et al 1996) for the 3-metric 7^^ = (fi'^fij, where 
(p is the conformal factor and fij the flat 3-metric in spherical coordinates. This 
approximation uses the maximal slicing condition and quasi-isotropic coordinates as 
gauge conditions. Under this approximation the resulting system consist in a hierarchy 
of elliptic equations (Cordero-Carrion et al 2009). 

To be able to study numerical dissipation effects we need a physical counterpart 
to calibrate our results. We use for this purpose the bulk viscosity of the fluid. The 
energy momentum tensor of a fluid with bulk viscosity (Ehlers 1961) is 

T^"" = p{l + e)u^'u'' + (P - Ce)/^^^ (2) 

where p is the rest-mass density, e is the speciflc internal energy, u'^ is the 4-velocity 
of the fluid, P is the pressure, C is the bulk viscosity, Q = ufj, is the expansion of 
the fluid and h^'^ = g^'^ + u^u'^ . The bulk viscosity appears as an isotropic term in 
the energy-momentum tensor in a very similar way to the pressure. Therefore we can 
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define a modified pressure P = P — QQ such that the energy-momentum tensor has 
the same form as a perfect fluid 



(3) 



where h = h — C^/p, being h = 1 + e + P/ p the rclativistic specific enthalpy. Using 
this form of the energy- momentum tensor is easy to modify an existing hydrodynamics 
code to include bulk viscosity, specially if the next assumptions are taken into account: 
i) we approximate the expansion assuming a post-Newtonian expansion and small 
perturbations, 

e^V-u + 0{c-^)+0{v^), (4) 

where ©(c^^) corresponds to first post-Newtonian corrections. All the simulations 
that we run with physical bulk viscosity in this work belong to this regimen, ii) The 
bulk viscosity coefficient ( is sufficiently small, such that P > 0. And iii) we neglect 
the contrilnition of the bulk viscosity in the computation of the sound speed needed 
for our numerical scheme. 

The hydrodynamics equations with bulk viscosity can be cast as a system of 
conservation laws (cf. Ibafiez et al 2000) 

(5) 
(6) 
(7) 



= 0, 
1 



where D = pW, Sj = phW^vj, r = phW^ — P — D are the conserved variables, 
v*^ = dx^/dt is the coordinate 3- velocity, = (w**— /3*)/a is the 3- velocity as measured 
by an Eulerian observer, and W = l/y'l — jijv^v^ the Lorentz factor. In the non- 
relativistic limit the eqiiations for a classical viscous ffuid can be recovered (Landau 
& Lifschitz 1987) which for constant ( result in the Navier-Stokes equations. 

We solve the coupled system of CFC spacetime evolution and GR hydrodynamics 
equations using the numerical code COCONUT (Dimmelmeier et al 2002, 2005). 
The numerical code uses standard high-resolution shock-capturing schemes for the 
hydrodynamics evolution in spherical polar coordinates, and spectral methods for the 
spacetime evolution. 



3. Estimating numerical viscosity 

To estimate the numerical viscosity of our code we have designed a simple spherical 
test. We consider a spherical fluid system of radius R and constant density in a static 
Minkowski spacetime. This system allow for discrete modes of radial oscillations (see 
appendix A) . Since the eigenfunctions of the modes are known it is possible to excite 
very accurately single modes of the system and follow their evolution. We evolve the 
system using a polytropic equation of state, therefore, since the system is adiabatic, 
there are no possible energy losses and the oscillations should keep constant amplitude 
as long as non-linear effects does not appear. Hence, any damping observed in the 
numerical simulation has to be caused by numerical dissipation effects. 

We have chosen a system with R = 1 and initial density po = 1. The equation 
of state is a polytrope of the form P = Kp^ with T = 4/3 and K = 1/3 x 10~^. 
We use a perturbation A{r) (see appendix A) of the velocity with an amplitude of 
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Figure 1. Time evolution of the first oscillation mode in the spherical test system. 
All evolutions shown are computed using the Marquina flux formula. Left panel 
shows the evolution of the perturbation a normalized to its initial value ao for a 
radial resolution rir = 80 and different reconstruction procedures: constant (black 
line), minmod (blue line) and MC (red line). The right panel shows the evolution 
of the amplitude of the perturbation Aa normalized to its initial value Aao using 
minmod reconstruction and different resolutions: 10 (black), 20 (blue), 40 (red), 
80 (green) and 160 (cyan). We also plot the fit of the lowest resolution simulation 
to an exponential decay (black line). 

v'q = 10~^ corresponding to the lowest frequency mode, uj = 0.094. This amplitude 
is sufficiently small for the oscillations to be considered linear. We have used second 
order Runge-Kutta method for the time evolution in all cases and two different flux 
formulae, Marquina (Donat el al 1998) and HLL (Harten & van Leer 1983). We have 
computed all models using different reconstruction techniques (constant and linear 
with minmod or MC slope limiters) and different resolutions (n^ = 10, 20, 40, 80 and 
160). We have also computed the models without physical bulk viscosity [C, — 0) and 
with non-zero values (C = 10"'', 10~^ and 10~^). We have evolved the system for 
2700 cj^^. To follow the evolution of the oscillations we computed the quantity 

a{t) = j dx^A{r)v''{r,t). (8) 

Left panel of figure [1] shows the time evolution of a/ao for different reconstruction 
procedures, oo being the initial value of a. The oscillation frequency coincides with the 
predicted by the linear analysis within 0.1% forrir = 80 . It can be seen that the order 
of the reconstruction has a strong influence in the numerical damping of the system, 
being much stronger for first order reconstruction (constant) than for second order 
(minmod or MC). In order to quantify this effect we estimate the damping time for 
each simulation as it is described next. In every case, we first compute the amplitude 
of a as the difference between two consecutive oscillations, Aa. In the right panel of 
figure [T] we plot the evolution of the amplitude for the minmod reconstruction case for 
five different resolutions. For very low resolution, the amplitude of the perturbation 
falls bellow the round-off error of the code within the duration of the simulation and 
it saturates. We fit next the the amplitude of the oscillations to an exponential decay 
Aa(t) = Aa(0) e"^*/"^ for each model. This procedure allow us to compute the value 
of the damping time of the amplitude, t/2. Since the energy of the perturbations 
scales quadratically with the velocity, the damping time of the energy is r. Finally, 
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Figure 2. Dependence of the numerical viscosity f with the resolution in the case 
of the test sphere. The left panel shows the case without physical bulk viscosity 
for different flux formulae, Marquina (filled circles) and HLL (open circles), and 
different reconstruction schemes, constant (back), minmod (red) and MC (green). 
The right panel shows the case with Marquina flux formula and MC reconstruction 
for different values of the physical viscosity C, = 10~* (red) 10~^ (green), 10~^ 
(blue) and (black). Dashed lines of the same color correspond to the value of 
the physical viscosity in each case. 



using (jB.9p it is possible to compute the physical bulk viscosity from the value of r as 
C = ^, (9) 

UJ T 

where (^yn = "^Pocl /uj and Cg is the sound speed. The left panel of figure [5] shows 
the behavior of the numerical viscosity with the radial resolution for the different 
numerical schemes tested in simulations without physical bulk viscosity (C = 0). To 
check the scaling with resolution we fit the viscosity values to a power law 

Cdyn ( A ) ' ^ ^ 

where A = 2TrCsUj^^ is the wavelength of the oscillation mode, being Cg the speed 
of sound. The results of the fitted parameters S and p for the different numerical 
methods are shown in table [T] In all cases S is of order one, although in general the 
Marquina flux formula provides about 30% less bulk viscosity than the same simulation 
with HLL. The scaling does not change substantially with the flux formula but only 
with the reconstruction scheme: for constant reconstruction we recover first order 
convergence and for linear reconstruction (minmod and MC) second order. These fits 
allow us to compute the numerical bulk viscosity of different numerical schemes. 

In the right panel of figure [5] we consider the simulations with physical bulk 
viscosity (C ^ 0). In this case the viscosity of the code decreases with resolution 
with a similar scaling as in the case without physical viscosity. For sufficiently high 
resolution, the viscosity converges to the value of the added physical viscosity. Our 
method to estimate the bulk viscosity overestimates, in the convergent regime, the real 
physical viscosity by ~ 25%. We have checked that the damping time and therefore 
the viscosity is not affected for a wide range of amplitudes of the perturbation (0.1- 
10~^) and hence we can conclude that non-linear effects does not play any role in 
the damping observed during this test. We have also checked that the results are 
independent of the CFL factor used for the time-step, by changing its value between 
0.8 and 10^"^. Since our final aim is to perform multidimensional simulations we have 
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Table 1. Coefficients S and p of the fitted numerical viscosity in radial oscillations 
of a test fluid sphere, for different flux formulae and reconstruction schemes. 



Flux formula reconstruc. p S 



Marq. 


constant 


0.97 


0.74 


HLL 


constant 


0.98 


1.06 


Marq. 


minmod 


2.01 


1.91 


HLL 


minmod 


2.03 


2.66 


Marq. 


MC 


2.23 


0.96 


HLL 


MC 


2.23 


1.28 



also performed some axisymmetric 2D models of this spherical system. We have chosen 
two representative radial resolutions Ur = 20 and 80 and varied the angular resolution 
ng = 8, 16 and 32. In all cases the results are indistinguishable from the ID spherical 
simulations with the same number of radial points. 

4. Numerical damping of neutron stars 

We apply our method to estimate viscosity in the case of radial oscillations of non- 
rotating neutron stars. Our initial model is a non-rotating rclativistic equilibrium 
configuration (Tolnian 1939, Oppcnhcimcr & Volkoff 1939). Wc use a polytropic 
equation of state P = K with adiabatic index F = 2 and polytropic constant 
K = 100 (units of G = c = Mq — 1). We choose central rest mass density to 
be pc = 7.9 X 10^** g cm~^ and gravitational mass Mg = \AMq. The resulting 
circumferential radius is Rc — 14.16 km, and the isotropic radial coordinate at the 
surface of the star is i? = 12.0 km. We evolve the system in dynamic spacetime 
adding an initial perturbation of the radial velocity corresponding to the fundamental 
radial oscillation mode with an amplitude of 10^'^. To compute the eigenfunction 
of the fundamental mode necessary for the perturbation, we use the eigenfunction 
recycling technique described in Dimmelmeier et al 2006. We perform ID spherical 
simulations with different radial resolutions of the fluid grid, = 80, 160, 320 and 
640, which is equally spaced from the center to 14.4 km. The time-step is computed 
using the CFL condition for the eigenvalues of the hydrodynamics system, resulting 
in = 1.13 X 10~^ ms for a grid with = 80 and a CFL factor 0.8. We use 
an artificial atmosphere to treat the vacuum surrounding the star as described in 
Dimmelmeier et al 2002. The treatment consists of resetting all numerical cells with 
p bellow a certain threshold to the value patm and setting the velocity to zero. The 
values for the threshold and the atmosphere are lO^*' and 10~^° times the initial 
central density respectively. We have checked that the results presented here do not 
change if we decrease the values of these two quantities. We use three radial domains 
for the spectral metric solver: two including the star and one compactified domain in 
the exterior. Each of the domains has either 17 or 33 collocation points for our low 
and high metric resolution simulations. Since the CFC approximation only contains 
elliptic equations, which are not restricted by the CFL condition, it is not necessary 
to solve the metric as frequently as the hydrodynamics equations. We use a metric 
computation rate of 71^/8, i.e. we compute the metric once every n^/S time step of 
the hydrodynamics. We use a parabolic extrapolation of all the metric quantities 
between consecutive metric computations, which has been shown to provide sufficient 
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Figure 3. Dependence of the numerical viscosity C, on the resolution in neutron 
star evolution. We show values for different flux formulae, Marquina (filled circles) 
and IfLL (open circles), and different reconstruction schemes, constant (back), 
minmod (red) and MC (green). In all cases the highest metric resolution case is 
plotted. The dashed line represents the function C/Cdyn = Ar/A. 



accuracy during the evolution (Dimmelmeier et al 2002). We do not add physical bulk 
viscosity in any of the simulations (C = 0), therefore the only dissipative processes in 
the evolution are of numerical nature. To compute the numerical bulk viscosity we 
use ([9]) as in the previous section but with ^dyn = 2pQ(?^/uj according to (|B.14p . where 
p is estimated as the central density. Since the numerical viscosity can depend on 
the location in the star the resulting value represents an average numerical viscosity. 
For our particular case and the fundamental radial mode (w = 9.11 x 10'^ Hz) the 
resulting value is Cdyn = 3.99 x 10'^^ g cm~^ s~^. This procedure gives us an idea of 
the order of magnitude of the averaged numerical viscosity on the star. More accurate 
computations of the damping-time could be done using the expressions of the appendix 
B or the procedure described in Cutler et al 1990. 

Figure [3] shows the variation of the viscosity with the radial resolution normalized 
to the wavelength of the perturbation computed at the center of the star. In general 
the viscosity decreases with increasing resolution. However the lowest resolution case 
with the minmod reconstruction shows an anomalously low numerical viscosity, which 
we think is an artifact of the low resolution. We plot the function C/Cdyn = Ar/A 
for comparison. The general trend is roughly first order decreasing the order for the 
highest resolution models, although we are using both first and second reconstruction 
schemes. The reason for the lower order of convergence is that all reconstruction 
schemes reduce to first order at discontinuities. Since we encounter a discontinuity at 
the surface any evolution will inevitably lead to first order convergent results. This 
is also a strong indication that the main numerical dissipation process in the neutron 
star evolution is probably due to the surface of the star, and therefore increasing the 
reconstruction order does not decreases the mean numerical viscosity responsible for 
the damping of global oscillations. 

We have checked the influence of the CFL factor, the metric calculation rate and 
the 2D effects for simulations with Ur = 80 radial points, Marquina fiux formula and 
MC reconstruction. One of the main differences with respect to the spherical test case 
of the previous section is that the equations that we are solving have sources due to the 
gravity terms. In general the presence of non-zero sources can increase the stiffness of 
the problem leading to inaccurate solutions, even unstable for extremely stiff sources. 
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We expect some influence of the CFL factor in our simulations since the reduction of 
the time step tends to cure the stiffness problems. We have performed simulations 
with CFL factor 0.8 (standard for all previous simulations), 0.1, 0.01 and 0.001, with a 
metric computation rate of 10 x 0.8/(CFL factor), which maintain constant the metric 
computation rate per unit evolution time. The amplitude of the oscillation at the 
end of the simulation varies about 0.5% at most. However most of the variation is 
observed between CFL factor 0.8 and 0.1, and decreasing the CFL factor further does 
not significantly changes the solution. Therefore the effect on the computed damping 
rate can be neglected. We find that the stiffness of the sources depends on how 
often the metric is computed. We also perform the same simulations with a metric 
computation rate 10. In this case the metric is computed more often per time unit 
than the previous one, as the CFL factor decreases. We find variations of about 25% 
in the final amplitude depending on the CFL factor used, although these variations 
depend strongly on the radial resolution. For = 160 the maximum variation due to 
CFL factor changes is smaller than 10%. Since the effect of the CFL factor becomes 
larger lowering the metric computation rate, we can conclude that the sources become 
more stiff if the metric is computed more often. An explanation for this effect is that if 
the metric is computed less often, the sources vary due to the space-time evolution on 
time-scales longer, than the hydrodynamic variables, which reduces the stiffness. We 
conclude that the stiffness of the sources is not significantly affecting the computation 
of the damping rates for the metric computation rate that we use in our regular 
simulations. However special care should be taken in simulations in which the metric 
computation rate is high (close to one), e.g. in black hole formation simulations, since 
the stiffness of the sources could lead to an increasing of the numerical viscosity of the 
code if the CFL factor is not appropriately lowered. We also perform 2D-axisymmetric 
simulations with 4, 8 16 and 32 angular grid points. Since the initial perturbation is 
radial, non-radial modes are not excited. We find an error in the angular velocity which 
is of order /v^ ~ 10~^, and which is closely related to the accuracy of numerical 
recovery of the primitive variables. Since the angular extension of the grid restricts 
the time step of the simulation in a factor proportional to l/ne, we find a similar 
behaviour increasing ue as it is found decreasing the CFL factor. 

5. Discussion 

We conclude that it is possible to estimate dissipation effects in numerical simulations 
which produce an effect similar to physical bulk viscosity. In a very simplified 
the test fluid sphere of constant density, the numerical viscosity scales with the order 
of the reconstruction. In the application of our method to neutron stars simulations 
we find that the behaviour of the numerical viscosity is more complicated since the 
gravity sources also contribute to the numerical dissipation processes, and the presence 
of the surface lowers the scaling of the viscosity to first order. 

To understand why the numerical inaccuracies can be modelled as a bulk viscosity 
term we consider the generic form of the flux formula in a Riemann solver 

F4^ = F' + AAi7, (11) 

where -Fnum is the numerical flux that is used in the time update, is an averaged 
value of the flux, A[/ represents the discontinuity at each cell interface and A the 
eigenvalue matrix. The particular form of the average and the matrix U depends on 
the particular Riemann solver. If a reconstruction scheme of order p is used, the value 
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of U at both sides of the interface in smooth regions of the flow agrees up to the 
p — 1-th derivative, therefore for any function f{U) 

Af{U)(x{ArfVPU (12) 

where represents p-th spatial derivatives. If one apphes this expression to the 
momentum equation 

dt iV^Sj) + d,F' + p {Ar)P V^P+^^S, = Q, (13) 

where depends on the specific Riemann solver used and Q represent the sources 
(r.h.s of (jA.2[) '). The new term appearing in the numerical version of the equations is 
a dissipative term since in includes ap+l derivative, while /i (Ar)f is the corresponding 
dissipative coefficient. For first order reconstruction, the dissipation term includes 
second derivatives that resemble to a viscosity term with viscosity proportional to Ar. 
For higher order reconstruction the new term resembles hyperviscosity coefficients. 
Note that with this interpretation we can identify the computed value of C in our 
simulations with a viscosity or hyperviscosity term, which scales with resolution with 
the order of the reconstruction scheme. Furthermore, we can argue that, since the new 
terms in (jl3p are responsible for both the bulk viscosity and shear viscosity terms, 
the value of the numerical shear viscosity of the code has to be of the same order 
of magnitude as the numerical bulk viscosity estimated in this work. If this case 
were true, the method described in this paper could be a powerful tool to estimate 
numerical dissipative effects in multidimensional simulations . However the method 
has some limitations due to the approximations that we considered. The integral 
expressions (|B.12p and (jB.lSp as well as the approximate expression for the expansion 
Q rely on the facts that a post-Newtonian expansion is possible and that velocity 
perturbations are small. This provides an order-of-magnitude estimate in the case 
of systems involving neutron stars or proto-neutron stars since 0{c~^) ^ 0.15. In 
this scenarios the velocity involved in oscillations is typically smaller than 0.1. In the 
vicinity of black holes, the post-Newtonian expansion is not convergent anymore, and 
the method can lead to large inaccuracies. Similar thing happens in the case of flows 
with a high Lorentz factor as those observed in jets. 

From our simulations we also conclude that the gravity can be an important 
source of numerical viscosity which has to be estimated in an appropriate way. The 
stiffening of the sources in the presence of rapidly changing spacetime, can lead to 
a strong increasing in the numerical viscosity if the CFL condition is not adapted 
accordingly. Although this effect is not a problem in the simulations presented in this 
work, it could be in case in which the spacetime evolves in similar time-scales as the 
fluid, e.g. in the formation of a black hole. It is therefore important to test in the 
future which are the real effects of numerical viscosity in such simulations. 

Finally, we note that all the viscosity results given in this paper are normalized 
to Cdyn, which dcpcuds on the frequency of the oscillation u. For a given astrophysical 
scenario, and the same numerical resolution, different modes will be affected in a 
different way depending on the frequency of the mode. Therefore, in order to estimate 
the resolution needed to evolve a system for a given amount of time with reasonably 
small damping one has to use the mode with higher frequency, of those who may be 
important in the dynamics of the system. Deciding which is this mode may be non 
trivial in non-linear simulations, since a strong numerical damping in high frequency 
modes, which are coupled to the lower frequency modes, can modify the damping 
times for the low frequency modes too. 
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Appendix A. Radial oscillations of an homogeneus sphere. 

We consider a spherically symmetric barotropic fluid in Minkowski spacetime, with 
equation of state P = P{p). The hydrodynamics equations (lA.ll) and (IA.2I) in this 
case read: 

dt {pW) + ^dr [r^pWv] = 0, (A.l) 

dt (^phW^v^ + ^dr r^phW^v^ = - drP, (A.2) 

where v = . The equilibrium solution of this system is trivially zero velocity, vq = 0, 
and constant pressure, P — Pq, and hence constant rest mass density, p — pq. We 
consider perturbations of the rest mass density and velocity, p — po + p' and v = v' , 
and hence 

P = Po+ ^p' = Po + hoclp'. (A.3) 
opo 

where is the speed of sound of the equilibrium model. The linearized hydrodynamics 
equations read 

dtv' + ^drp' - ^drO - 0, (A.4) 
Po Poho 

dtp' + podrv' + -pov' - 0, (A.5) 
r 

If ( is sufficiently small, i.e. C << ^oPOj the viscosity terms can be neglected in the 
computation of the perturbations spectrum. In this case equations (|A.4p and (jA.Sp 
can be combined in a wave equation for p' 

dttp' - c%rP' - —drp' = 0. (A.6) 
r 

The system of equations (IA.4I) and (|A.5p admits oscillatory solutions of the form 
v' = A{r) smuit and p' — B{r) cosujt and hence (|A.6p results in 

2 uj"^ 

drrB + -drB + B = 0. (A.7) 
r cj 

This equation has the form of the Lane-Emden equation of index 1, which have 
solutions regular at r = (Chandrasekhar 1967) of the form 

^ ; sin (kr) , . ^ 

B{r)=p'^^-1, (A.8) 

where k = ui/cs and is a parameter which controls the amplitude of the density 
perturbation. From (|A.4p the solution for the velocity perturbation is 

. , , . kr cos (k r) — sin (kr\ ... 
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where Vq controls the amphtude of the velocity perturbation and is related to p'^ by 

^ = (A.IO) 

If we impose boundary conditions at the surface v'{R) = 0, i.e., A{R) = we find a 
discrete spectrum of modes given by the solutions of 

tan (fc i?) = fc i?. (A.ll) 

We have computed the roots of this equations by means of a bisection algorithm. The 
first five numerical solutions correspond to k R — 4.49, 7.73, 10.90, 14.07 and 17.11. 

Appendix B. Damping time of an oscillating spherical star 

The rate of change of the energy due to bulk viscosity of a pulsating star is (cf. Cutler 
et al 1990) 

'^ = ^An j\rr^cj,^(:\Q\\ (B.l) 
and, provided that the energy of the pulsations E is known, the damping time is 



r--2El^-J . (B.2) 

where <> denotes the time average over a cycle. 

The energy stored in the radial oscillations of a spherical star can be computed 
as the energy difference between the star in equilibrium and the perturbed system. 
Previous computations (Meltzer & Thorne 1966, Glass & Lindblom 1983) used 
Schwarzschild coordinates in their computations. Instead of performing the coordinate 
transformation to the choice of the present work we find it easier and more instructive 
to compute the energy of the oscillations directly in our coordinates, although the 
result should be identical. 

The ADM energy is a conserved quantity which in spherical symmetry and under 
our gauge choice is 

Eabm = -2 J dr r^A^ = in J dr r^cb'' iphW^ -P+ -j^j , (B.3) 

being A the Laplacian with respect to the flat 3-metric and Kij the extrinsic curvature 
of the induced 3-metric 7^ . Since the extrinsic curvature vanishes for the equilibrium 
system in our gauge choice the equilibrium ADM energy is 

Eadmo = 4^ / drr^<j>l{poh„ - Pq)- (B.4) 
Jo 

The ADM energy does not change with time and hence we can compute its value 
by evaluating the integral at the oscillation phase with maximum velocity. In this 
phase, due to the continuity equation (IA.1[) . the variation of ip^pW with respect to 
the equilibrium is zero. The leading term in the perturbation corresponds to quadratic 
terms in the velocity. The energy of the oscillations is thus E = Eabm — Eabm which 
results in 

E = 47: I dr r^^s Poho^v'' - {poho - 5Po)^ + " (B.5) 



/o 



Numerical viscosity in GR hydrodynamics 



13 



for adiabatic perturbations. Note that here 4>' = (j) — 4>q corresponds to the variation 
of (f> with respect to the equihbrium for the phase in which v' is maximum and hence 
is a term quadratic in v' . 

If we apply this expressions to the case of the spherical fluid in Minkowski 
spacetime of the appendix A the resulting expressions are 

E^^^j^ drr^ipo^W' =Pof'o^^sin2(fci?), (B.6) 

= -Cv'lirRsin^ (kR), (B.7) 

where we have explicitly used that (jA.ll[) is fulfilled to perform the integration and 
that < |ep >= (VA)2/2 = uj^B^/{2pl). The damping time is therefore 

(B.8) 

It is convenient to express it as 

Sdvn 

where Cdyn = 2/9oCs/w. For C ^ Cdyn the damping of the mode occurs in dynamical 
time-scales, while if C << Cdyn the damping occurs in secular timescales. 

In the case of a perturbed relativistic star the computation of the damping rate 
can not be computed analytically. However we can still make an order of magnitude 
estimation. For this purpose it is useful to truncate the equations (jB.ip and (IB. 51) to 
the leading order in the post-Newtonian expansion 

E =4tt drr^-pov'^{l + 0{c-^)), (B.IO) 

.In 2 



(B.9) 



dE 



R 



^=-47r/ drr^Ciep (1 + 0(0-2)) . (B.ll) 
dt Jo 

which corresponds to the Newtonian expression for the kinetic energy, beging 0{c~^) ~ 
0{v'^) ~ 0{GM / R) the first post-Newtonian corrections (see e.g. Blanchet et al 1980). 
If we assume dependence of the perturbations, where k is the wavenumber, 

then < |6p >~ k'^ < v'^ >— k'^v'^^/2. The energy and energy losses result in this 
case 

fR 

E = 27: drr2po<ax (l + (^(c"')) , (B.12) 
Jo 

^ -2nk^cj\rr\l^^ (l + 0{c-')) , (B.13) 



and hence 



2a;po _ 2poCs _ Cdyn 



(B.14) 



k\ Cu: ' C ' 
where p is an averaged density weighed by the eigenfunction A{r) and Cdyn = '^Po(?sl^- 
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